
dir.create("HRS_eventstudy")

#already excluded women, highly educated, ever self-employed, ever on SSI alone, or who never work while healthy
hrsdata<-as.data.table(read.dta13(paste0("hrs_health_",event,".dta")))


hrsdata[,homeowner:=as.numeric(homevalue==0)]
hrsdata[,veteran:=as.numeric(veteran)-1]
hrsdata[,baseyear:=-2]
hrsdata[is.na(basediffyear),basediffyear:=2]
setnames(hrsdata,c("cohort","pweight"),c("birthcohort","personweight"))
hrsdata[,work:=as.numeric(work)-1]
hrsdata[,sp_work:=as.numeric(sp_work)-1]
hrsdata[,race:=as.numeric(race)]
hrsdata[,birthcohort:=as.numeric(birthcohort)]
#hrsdata<-hrsdata[!is.na(race),]
hrsdata[,empyet:=cumsum(work)>=1,by=id]
hrsdata[,initwork:=work]
hrsdata[,diffcheck := basediffyear==2]
hrsdata[,basecheck:=diffcheck==1 & empyet == 1 & 
          (year <= firstrecyear | is.na(firstrecyear)) & 
          (year <= firstappyear | is.na(firstappyear))]

hrsdata[is.na(anyonset_age),anyonset_age:=Inf]
hrsdata[is.na(refonset_age),refonset_age:=Inf]


hrsdata[,refonset_time:=year-refonset_year]
hrsdata[,prepartner:=haspartner[year==refonset_year],by=.(id)]
hrsdata[,sp_prework:=max(sp_work[refonset_time<0]),by=.(id)]
hrsdata[,sp_prewage:=max(sp_wage[refonset_time<0 & sp_prework==1]),by=.(id)]
sp_prewagevals <- hrsdata[prepartner==1 & sp_wage>0 & sp_work==1 & refonset_time < 0,max(sp_wage,na.rm=TRUE),by=.(id)]$V1
sp_prewagevals<-quantile(sp_prewagevals[!is.na(sp_prewagevals)])
hrsdata[,sp_prewagerank:=max(which(max(sp_prewage) > sp_prewagevals)),by=.(id)]
hrsdata[is.infinite(sp_prewagerank),sp_prewagerank:=NA]
hrsdata[is.na(prepartner),sp_prewagerank:=NA]
hrsdata[,sp_prewagerank:=sp_prewagerank+1]
hrsdata[prepartner==0,sp_prewagerank:=0]
hrsdata[prepartner==1 & sp_prework==0,sp_prewagerank:=1]

hrsdata[!is.na(sp_wage),sp_cumwage:=cummax(sp_wage),by=id]
hrsdata[!is.na(sp_work),sp_cumwork:=cummax(sp_work),by=id]
hrsdata[,row:=.I]
hrsdata[,sp_wagerank:=max(which(sp_cumwage > sp_prewagevals)),by=row]
hrsdata[sp_wagerank==5,sp_wagerank:=4]
hrsdata[is.infinite(sp_wagerank),sp_wagerank:=NA]
hrsdata[is.na(haspartner),sp_wagerank:=NA]
hrsdata[,sp_wagerank:=sp_wagerank+1]
hrsdata[haspartner==0,sp_wagerank:=0]
hrsdata[haspartner==1 & sp_cumwork==0,sp_wagerank:=1]

#Check: patterns in take-up by spousal wages consistent with story? Somewhat.
hrsdata[refonset_time>=0 & refonset_time <= 6,mean(onDI_rec),by=sp_prewagerank]
hrsdata[refonset_time>=0 & refonset_time <= 6,mean(yetapplyDI_rec),by=sp_prewagerank]

#employment (cumulative) and marital status of next year:
hrsdata[,empyet1:=c(empyet[-1],NA),by=id]
hrsdata[is.na(empyet1),empyet1:=empyet]
hrsdata[,empyet:=empyet1]
hrsdata[,empyet1:=NULL]
hrsdata[,marnow1:=c(haspartner[-1],NA),by=id]
hrsdata[is.na(marnow1),marnow1:=haspartner]
hrsdata[,haspartner:=marnow1]
hrsdata[,marnow1:=NULL]
hrsdata[,sp_wagerank:=c(sp_wagerank[-1],NA),by=id]
hrsdata[,N:=1]

wealthcutoffs<-hrsdata[age==refonset_age+baseyear & !is.na(wealth_nonh_real),quantile(wealth_nonh_real,probs=c(0.25,0.5,0.75))]
hrsdata[!is.na(wealth_nonh_real)&wealth_nonh_real<wealthcutoffs[1],rich:=0]
for(cut in 1:length(wealthcutoffs)){
  hrsdata[!is.na(wealth_nonh_real)&wealth_nonh_real>=wealthcutoffs[cut],rich:=cut]
}

covariates<-1
timevar = "age"
cohortvar = "refonset_age"
anycohortvar = "anyonset_age"

#For comparisons of means, restricting to valid treated units.
#If doing event study, I can't drop potential control units:
if(eventmeans==TRUE){
#already excluded women, highly educated, ever self-employed, ever on SSI alone, or who never work while healthy
#Just need to exclude to severe events where I can observe timing
hrsdata<-hrsdata[get(cohortvar)==get(anycohortvar),]

hrsdata[,time:=get(timevar) - get(cohortvar)]
hrsdata[,group:=as.factor(get(gvar)[time==baseyear]),by=.(id)]

hrsdata[,basecheck:=basecheck[time==baseyear],by=.(id)]
hrsdata<-hrsdata[basecheck==1,]

hrsdata[,agebin:=floor(age/10)*10]
hrsdata[,agebin:=max(agebin * (time==0)),by=id]
hrsdata<-hrsdata[!is.infinite(anyonset_age),]
}

